library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readr)
library(tf)
##
## Attaching package: 'tf'
##
## The following objects are masked from 'package:stats':
##
## sd, var
library(tidyfun)
library(ggplot2)
library(refund)
library(mgcViz)
## Loading required package: mgcv
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
## Loading required package: qgam
## Registered S3 method overwritten by 'mgcViz':
## method from
## +.gg ggplot2
##
## Attaching package: 'mgcViz'
##
## The following objects are masked from 'package:stats':
##
## qqline, qqnorm, qqplot
library(rgl)
library(plotly)
##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
library(fdaoutlier)
library(purrr)
library(patchwork)
library(table1)
##
## Attaching package: 'table1'
##
## The following objects are masked from 'package:base':
##
## units, units<-
sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Linux Mint 22.1
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/Berlin
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] table1_1.4.3 patchwork_1.3.0 fdaoutlier_0.2.1 plotly_4.10.4
## [5] rgl_1.3.18 mgcViz_0.2.0 qgam_2.0.0 mgcv_1.9-3
## [9] nlme_3.1-168 refund_0.1-37 tidyfun_0.0.98 tf_0.3.5
## [13] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
## [17] purrr_1.0.4 readr_2.1.5 tidyr_1.3.1 tibble_3.3.0
## [21] ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.4 bitops_1.0-9 gridExtra_2.3
## [4] rlang_1.1.6 magrittr_2.0.3 matrixStats_1.5.0
## [7] compiler_4.5.1 vctrs_0.6.5 pkgconfig_2.0.3
## [10] fastmap_1.2.0 backports_1.5.0 magic_1.6-1
## [13] promises_1.3.2 deSolve_1.40 rmarkdown_2.29
## [16] tzdb_0.5.0 pracma_2.4.4 nloptr_2.2.1
## [19] xfun_0.52 cachem_1.1.0 jsonlite_2.0.0
## [22] later_1.4.2 parallel_4.5.1 cluster_2.1.8.1
## [25] R6_2.6.1 bslib_0.9.0 stringi_1.8.7
## [28] RColorBrewer_1.1-3 GGally_2.3.0 extrafontdb_1.0
## [31] boot_1.3-31 rainbow_3.8 jquerylib_0.1.4
## [34] Rcpp_1.0.14 iterators_1.0.14 knitr_1.50
## [37] zoo_1.8-14 base64enc_0.1-3 hdrcde_3.4
## [40] extrafont_0.19 httpuv_1.6.16 Matrix_1.7-3
## [43] splines_4.5.1 timechange_0.3.0 tidyselect_1.2.1
## [46] rstudioapi_0.17.1 abind_1.4-8 yaml_2.3.10
## [49] viridis_0.6.5 doParallel_1.0.17 codetools_0.2-20
## [52] miniUI_0.1.1.1 RLRsim_3.1-8 lattice_0.22-5
## [55] plyr_1.8.9 shiny_1.10.0 ks_1.14.3
## [58] withr_3.0.2 S7_0.2.0 pbs_1.1
## [61] evaluate_1.0.3 ggstats_0.9.0 grpreg_3.5.0
## [64] mclust_6.1.1 pillar_1.10.2 fda_6.2.0
## [67] KernSmooth_2.23-26 checkmate_2.3.2 foreach_1.5.2
## [70] reformulas_0.4.1 pcaPP_2.0-5 generics_0.1.3
## [73] RCurl_1.98-1.17 hms_1.1.3 scales_1.4.0.9000
## [76] minqa_1.2.8 xtable_1.8-4 gamm4_0.2-7
## [79] glue_1.8.0 lazyeval_0.2.2 tools_4.5.1
## [82] fds_1.8 data.table_1.17.0 lme4_1.1-37
## [85] mvtnorm_1.3-3 grid_4.5.1 Rttf2pt1_1.3.12
## [88] rbibutils_2.3 colorspace_2.1-1 Formula_1.2-5
## [91] cli_3.6.5 viridisLite_0.4.2 gtable_0.3.6
## [94] sass_0.4.10 digest_0.6.37 htmlwidgets_1.6.4
## [97] farver_2.1.2 htmltools_0.5.8.1 lifecycle_1.0.4
## [100] httr_1.4.7 mime_0.13 MASS_7.3-65
strain_df<-readRDS("papillary muscle_data.rds")
strain_df$id<-factor(strain_df$id)
scalar_covariates<-readRDS("scalar_covariates.rds")
scalar_covariates$id<-factor(scalar_covariates$id)
# put irregular functions x that vary in length/start/end on a new_domain
# on a regular grid of length n_arg
tf_rectify <- function(x, new_domain = c(0, 1), n_arg = 30, ...) {
x_val <- tf_evaluations(x)
x_arg <- tf_arg(x) |> ensure_list()
# create new grids with same proportional spacings as old grids on the new domain:
new_arg_irreg <- purrr::map(x_arg,
function(x) (((x - min(x))/diff(range(x))) * # scale to [0,1]
diff(new_domain)) + # scale to [0, new_length]
new_domain[1]) # push to new domain
new_arg_grid <- seq(new_domain[1], new_domain[2], length = n_arg)
# put data on new common domain and re-eval on new grid
tfd(x_val, new_arg_irreg, ...) |> tfd(arg = new_arg_grid, ...)
}
# get regular evenly spaced time points for the curves
strain_df_clean <- strain_df %>%
mutate(
strain_str = ifelse(strain_curve %in% c("GLS4", "GLS2", "GLS3"), "GLS",
ifelse(strain_curve %in% c("LA4", "LA2"), "LA", strain_curve)
),
strain_reg = tfd(strain_val[, seq(0, 1, length.out = 60), interpolate = TRUE]),
id = factor(id)
)
# show complete dataset before aggregation:
ggplot(strain_df_clean) +
geom_spaghetti(aes(y = strain_reg, color = strain_str, linetype = cycle_number)) +
facet_wrap(~id) + labs(title = "Full Dataset")
# global mean strain curves
strain_df_clean |>
group_by(strain_str) %>%
summarize(
strain_mean = mean(strain_reg),
strain_se = sd(strain_reg) / sqrt(nlevels(id))
) %>%
filter(strain_str %in% c("APM", "PPM", "TOR", "LA", "GLS")) %>%
ggplot() +
geom_spaghetti(aes(y = strain_mean, color = strain_str)) +
geom_errorband(aes(
ymin = strain_mean - 2 * strain_se,
ymax = strain_mean + 2 * strain_se,
fill = strain_str
)) + labs(title = "Global Mean Strain Curves")
# average curves per patient: gls=mean(gls2,gls3, gls4), LA=mean (LA2, LA4)
strain_df_patientwise <- strain_df_clean %>%
group_by(id, strain_str) %>%
summarise(
strain_m = mean(strain_reg)
)
## `summarise()` has grouped output by 'id'. You can override using the `.groups`
## argument.
ggplot(strain_df_patientwise) +
geom_spaghetti(aes(y = strain_m, color = strain_str)) +
facet_wrap(~id) + labs(title = "Patient-wise Mean Strain Curves")
# putting in right format for analysis and adding the scalar covariates of interest
strain_df_patientwise_long <- strain_df_patientwise %>%
dplyr::select(c("id", "strain_str", "strain_m")) %>%
pivot_wider(names_from = strain_str, values_from = strain_m)
strain_df_patientwise_long <- strain_df_patientwise_long[-43, ]
strain_df_patientwise_long <- left_join(strain_df_patientwise_long, scalar_covariates, by = "id") #
strain_analysisdata <- strain_df_patientwise_long
strain_analysisdata_summaries <- strain_analysisdata %>%
mutate(
tor_max_time = TOR|>tf_where(value == max(value), "first"),
tor_max = TOR|>tf_fmax(),
apm_min_time = APM|>tf_where(value == min(value), "first"),
apm_min = APM|>tf_fmin(),
ppm_min_time=PPM|>tf_where(value == min(value), "first"),
ppm_min=PPM|>tf_fmin(),
gls_min_time=GLS|>tf_where(value == min(value), "first"),
gls_min=GLS |> tf_fmin())
table1(~
Age + factor(Gender) + factor(HTN) + factor(DIABETES) + weight + height +
heart_rate + tor_max + tor_max_time + gls_min + gls_min_time +
apm_min + apm_min_time + ppm_min + ppm_min_time,
data = strain_analysisdata_summaries)
| Overall (N=75) |
|
|---|---|
| Age | |
| Mean (SD) | 37.4 (14.9) |
| Median [Min, Max] | 36.2 [15.6, 78.1] |
| factor(Gender) | |
| 0 | 16 (21.3%) |
| 1 | 59 (78.7%) |
| factor(HTN) | |
| 0 | 63 (84.0%) |
| 1 | 12 (16.0%) |
| factor(DIABETES) | |
| 0 | 68 (90.7%) |
| 1 | 7 (9.3%) |
| weight | |
| Mean (SD) | 75.0 (13.3) |
| Median [Min, Max] | 76.5 [47.0, 115] |
| Missing | 5 (6.7%) |
| height | |
| Mean (SD) | 173 (8.21) |
| Median [Min, Max] | 175 [151, 191] |
| Missing | 5 (6.7%) |
| heart_rate | |
| Mean (SD) | 75.3 (13.7) |
| Median [Min, Max] | 72.5 [43.0, 128] |
| Missing | 1 (1.3%) |
| tor_max | |
| Mean (SD) | 25.2 (10.6) |
| Median [Min, Max] | 24.0 [9.57, 54.1] |
| tor_max_time | |
| Mean (SD) | 0.412 (0.0649) |
| Median [Min, Max] | 0.407 [0.271, 0.593] |
| gls_min | |
| Mean (SD) | -17.3 (2.20) |
| Median [Min, Max] | -17.2 [-22.8, -12.2] |
| Missing | 3 (4.0%) |
| gls_min_time | |
| Mean (SD) | 0.429 (0.0557) |
| Median [Min, Max] | 0.441 [0.288, 0.593] |
| Missing | 3 (4.0%) |
| apm_min | |
| Mean (SD) | -17.1 (3.80) |
| Median [Min, Max] | -17.0 [-29.7, -7.22] |
| apm_min_time | |
| Mean (SD) | 0.509 (0.0765) |
| Median [Min, Max] | 0.492 [0.322, 0.712] |
| ppm_min | |
| Mean (SD) | -17.8 (3.73) |
| Median [Min, Max] | -17.9 [-25.4, -8.63] |
| ppm_min_time | |
| Mean (SD) | 0.509 (0.0791) |
| Median [Min, Max] | 0.508 [0.322, 0.695] |
m3_data <- strain_analysisdata %>%
select(id, APM, PPM, TOR, Age, Gender, heart_rate)
# plotting the apm, ppm and torsion curves
f_summary <- m3_data %>%
dplyr::select("id", "APM", "PPM", "TOR") %>%
pivot_longer(!id, names_to = "Curve", values_to = "strain") %>%
group_by(Curve) %>%
summarise(
strain_mean = mean(strain),
strain_sd = sd(strain)
) %>%
mutate(
maximum = strain_mean |> tf_where(value == max(value), "first"),
minimum = strain_mean |> tf_where(value == min(value), "first")
)
f_summary %>%
ggplot() +
geom_spaghetti(aes(y = strain_mean, color = Curve)) +
geom_errorband(aes(
ymax = strain_mean + 2 * strain_sd / sqrt(nrow(m3_data)),
ymin = strain_mean - 2 * strain_sd / sqrt(nrow(m3_data)),
fill = Curve
)) +
geom_vline(xintercept = 0.41, color = "blue", linetype = 2) + # time of maximum torsion
geom_vline(xintercept = 0.50, linetype = 2) + # time of minimum PM (ie.maximum contraction)
labs(
x = "Cardiac cycle proportion",
y = "",
caption = "mean+/-2 standard errors",
legend = "Curve"
) +
theme_bw()
# comparing time of maximum PM contraction to time of maximum torsion occurrence
APM_min<-m3_data$APM |> tf_where(value == min(value), "first")
PPM_min<- m3_data$PPM |> tf_where(value == min(value), "first")
Tor_max<-m3_data$TOR |> tf_where(value == max(value), "first")
mean(c(APM_min,PPM_min))
## [1] 0.509
mean(APM_min)
## [1] 0.509
sd(APM_min)
## [1] 0.0765
mean(PPM_min)
## [1] 0.509
sd(PPM_min)
## [1] 0.0791
mean(Tor_max)
## [1] 0.412
sd(Tor_max)
## [1] 0.0649
(res1 <- t.test(APM_min, Tor_max, paired = TRUE))
##
## Paired t-test
##
## data: APM_min and Tor_max
## t = 11, df = 74, p-value <2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 0.0798 0.1146
## sample estimates:
## mean difference
## 0.0972
(res2 <- t.test(PPM_min, Tor_max, paired = TRUE))
##
## Paired t-test
##
## data: PPM_min and Tor_max
## t = 14, df = 74, p-value <2e-16
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## 0.0828 0.1111
## sample estimates:
## mean difference
## 0.0969
(res3 <- t.test(PPM_min, APM_min, paired = TRUE))
##
## Paired t-test
##
## data: PPM_min and APM_min
## t = -0.02, df = 74, p-value = 1
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.0183 0.0179
## sample estimates:
## mean difference
## -0.000226
# min apm/ppm ( maximum of papillary muscles contraction) occurs significantly later than maximum torsion
# association across patients:
layout(t(1:2))
plot(APM_min, Tor_max, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,1), col = "gold")
plot(PPM_min, Tor_max, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,1), col = "gold")
layout(t(1:2))
plot(APM_min, Tor_max - APM_min, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,0), col = "gold")
plot(PPM_min, Tor_max - PPM_min, asp = 1, col = rgb(0,0,0,.3)); abline(c(0,0), col = "gold")
# patients with later APM_min (PPM_min as well to lesser extent) tend to have larger delay between Tor_max and APM_min
# selecting only diastolic curve value, diastole begins at the maximum of torsion
m4_data <- strain_analysisdata %>%
select(id, TOR, APM, PPM, Gender, Age, heart_rate) %>%
mutate(tor_max = TOR |> tf_where(value == max(value), "first")) %>%
mutate(
tor_diastole = TOR |> tf_zoom(tor_max, 1),
apm_diastole = APM |> tf_zoom(tor_max, 1),
ppm_diastole = PPM |> tf_zoom(tor_max, 1)
)
## Warning: There were 222 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `tor_diastole = tf_zoom(TOR, tor_max, 1)`.
## Caused by warning:
## ! Combining incompatible <tfd_reg> with <tfd_reg> by casting to <tfd_irreg>.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 221 remaining warnings.
m4_data$heart_rate[is.na(m4_data$heart_rate)] <-
mean(m4_data$heart_rate, na.rm = TRUE) # imputing the mean for one missing heart rate observation
# obtaining regular and same domain over diastole torsion curves
m4_data <- mutate(m4_data,
tor_diast = tf_rectify(tor_diastole),
apm_diast = tf_rectify(apm_diastole),
ppm_diast = tf_rectify(ppm_diastole)
)
m4_data %>% ggplot() +
geom_spaghetti(aes(y = tor_diast)) +
labs(title = "Rectified diastole curves") +
m4_data %>% ggplot() +
geom_spaghetti(aes(y = apm_diast)) +
m4_data %>% ggplot() +
geom_spaghetti(aes(y = ppm_diast))
# data prep for pffr
m4_data$tor_d_mat <- m4_data$tor_diast |> as.matrix()
m4_data$apm_d <- with(m4_data, apm_diast - mean(apm_diast)) |> as.matrix()
m4_data$ppm_d <- with(m4_data, ppm_diast - mean(ppm_diast)) |> as.matrix()
# diastolic concurrent fof regression
m4 <- pffr(tor_d_mat ~ s(apm_d) + s(ppm_d) + Gender + Age + heart_rate,
yind = seq(0, 1, length.out = 30),
family = "gaulss", data = m4_data
)
summary(m4)
##
## Family: gaulss
## Link function: identity logb
##
## Formula:
## tor_d_mat ~ s(apm_d) + s(ppm_d) + Gender + Age + heart_rate
##
## Constant coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.783 0.729 1.07 0.28
## (Intercept).1 1.016 0.024 42.34 <2e-16
##
## Smooth terms & functional coefficients:
## edf Ref.df Chi.sq p-value
## Intercept(yindex) 10.21 19.00 56.0 < 2e-16
## s(apm_d) 19.03 22.73 276.4 < 2e-16
## s(ppm_d) 13.48 17.23 60.2 2.4e-06
## Gender(yindex) 1.01 1.02 0.0 1
## Age(yindex) 3.18 3.54 49.9 < 2e-16
## heart_rate(yindex) 3.69 3.97 98.8 < 2e-16
## NA 18.80 19.00 2264.5 < 2e-16
##
## R-sq.(adj) = Deviance explained = 100%
## -REML score = 5896.4 Scale est. = 1 n = 2250 (75 x 30)
pffr.check(m4, rep = 100)
##
## Method: REML Optimizer: outer newton
## full convergence after 14 iterations.
## Gradient range [-0.0017,0.000295]
## (score 5896 & scale 1).
## Hessian positive definite, eigenvalue range [0.00169,6.2].
## Model rank = 125 / 125
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(yindex.vec) 19.00 10.21 0.86 <2e-16
## ti(apm_d,yindex.vec) 35.00 19.03 1.07 1
## ti(ppm_d,yindex.vec) 35.00 13.48 1.07 1
## s(yindex.vec):Gender 5.00 1.01 0.86 <2e-16
## s(yindex.vec):Age 5.00 3.18 0.86 <2e-16
## s(yindex.vec):heart_rate 5.00 3.69 0.86 <2e-16
## s.1(yindex.vec) 19.00 18.80 0.86 <2e-16
m4_gam <- m4
class(m4_gam) <- class(m4)[-1]
e <- getViz(m4_gam)
qq.gamViz(e, level = .9, CI = "quantile")
matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T) |>
cov() |>
filled.contour()
plotly::plot_ly(
z = ~ cov(matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T)),
type = "surface",
colorscale = "Viridis", # Color scheme - you can change this
showscale = TRUE # Show color scale legend
)
matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T) |>
cor() |>
filled.contour()
plotly::plot_ly(
z = ~ cor(matrix(residuals(m4, reformat = FALSE), nrow = m4$pffr$nobs, byrow = T)),
type = "surface",
colorscale = "Viridis", # Color scheme - you can change this
showscale = TRUE # Show color scale legend
)
# changed color scheme here so we can easier distinguish between
# time-value regions with no significant association to TOR (white) and
# small-effect-but-still-statistically-significant regions (grey-ish hues)
plot(sm(e, 2)) +
l_fitRaster(pTrans = function(.p) .p < 0.05) +
scale_fill_gradient2(low = "blue", mid = "lightgrey", high = "red") +
l_fitContour(bins = 20) + labs(y = "Diastole proportion", x = "APM", title = "") +
guides(fill = guide_legend(title = "Effect on untwist"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
plot(m4, scheme = 1, select = 2, ticktype = "detailed")
plot(sm(e, 3)) +
l_fitRaster(pTrans = function(.p) .p < 0.05) +
scale_fill_gradient2(low = "blue", mid = "lightgrey", high = "red") +
l_fitContour(bins = 20) + labs(y = "Diastole proportion", x = "PPM", title = "") +
guides(fill = guide_legend(title = "Effect on untwist"))
## Scale for fill is already present.
## Adding another scale for fill, which will replace the existing scale.
plot(m4, scheme = 1, select = 3, ticktype = "detailed")
diastole_summary <- m4_data %>%
pivot_longer(cols = 12:14, names_to = "Curve", values_to = "strain") %>%
group_by(Curve) %>%
summarise(
strain_mean = mean(strain),
strain_sd = sd(strain)
)
diastole_summary %>% ggplot() +
geom_spaghetti(aes(y = strain_mean, color = Curve)) +
geom_errorband(aes(
ymax = strain_mean + 2 * strain_sd / sqrt(nrow(m4_data)),
ymin = strain_mean - 2 * strain_sd / sqrt(nrow(m4_data)),
fill = Curve
)) +
labs(
x = "Diastole proportion",
y = "",
caption = "mean+/-2 standard errors",
legend = "Curve"
) +
theme_bw()